home *** CD-ROM | disk | FTP | other *** search
/ The Programmer Disk / The Programmer Disk (Microforum).iso / xpro / c2 / pro3 / fftr.c < prev    next >
Text File  |  1987-06-20  |  4KB  |  235 lines
  1. /*                            fftr.c
  2.  *
  3.  *    FFT of Real Valued Sequence
  4.  *
  5.  *
  6.  *
  7.  * SYNOPSIS:
  8.  *
  9.  * double x[], sine[];
  10.  * int m;
  11.  *
  12.  * fftr( x, m, sine );
  13.  *
  14.  *
  15.  *
  16.  * DESCRIPTION:
  17.  *
  18.  * Computes the (complex valued) discrete Fourier transform of
  19.  * the real valued sequence x[].  The input sequence x[] contains
  20.  * n = 2**m samples.  The program fills array sine[k] with
  21.  * n/4 + 1 values of sin( 2 PI k / n ).
  22.  *
  23.  * Data format for complex valued output is real part followed
  24.  * by imaginary part.  The output is developed in the input
  25.  * array x[].
  26.  *
  27.  * The algorithm takes advantage of the fact that the FFT of an
  28.  * n point real sequence can be obtained from an n/2 point
  29.  * complex FFT.
  30.  *
  31.  * A radix 2 FFT algorithm is used.
  32.  *
  33.  * Execution time on an LSI-11/23 with floating point chip
  34.  * is 1.0 sec for n = 256.
  35.  *
  36.  *
  37.  * PROGRAM AUTHOR:
  38.  *
  39.  * S. L. Moshier
  40.  *
  41.  *
  42.  * REFERENCE:
  43.  *
  44.  * E. Oran Brigham, The Fast Fourier Transform;
  45.  * Prentice-Hall, Inc., 1974
  46.  *
  47.  */
  48. /*            fftr.c            */
  49.  
  50.  
  51. static short n0 = 0;
  52. static short n4 = 0;
  53. static short msav = 0;
  54.  
  55. extern double PI;
  56.  
  57.  
  58. fftr( x, m0, sine )
  59. double x[];
  60. int m0;
  61. double sine[];
  62. {
  63. int th, nd, pth, nj, dth, m;
  64. int n, n1, n2, j, k, l, r;
  65. double xr, xi, tr, ti, co, si;
  66. double a, b, c, d, bc, cs, bs, cc;
  67. double *p, *q;
  68. double sin();
  69.  
  70. /* Array x assumed filled with real-valued data */
  71. /* m0 = log2(n0)            */
  72. /* n0 is the number of real data samples */
  73.  
  74. if( m0 != msav )
  75.     {
  76.     msav = m0;
  77.  
  78.     /* Find n0 = 2**m0    */
  79.     n0 = 1;
  80.     for( j=0; j<m0; j++ )
  81.         n0 <<= 1;
  82.  
  83.     n4 = n0 >> 2;
  84.  
  85.     /* Calculate array of sines */
  86.     xr = 2.0 * PI / n0;
  87.     for( j=0; j<=n4; j++ )
  88.         sine[j] = sin( j * xr );
  89.     }
  90.  
  91. n = n0 >> 1;    /* doing half length transform */
  92. m = m0 - 1;
  93.  
  94.  
  95. /*                            fftr.c    */
  96.  
  97. /*  Complex Fourier Transform of n Complex Data Points */
  98.  
  99. /*    First, bit reverse the input data    */
  100.  
  101. for( k=0; k<n; k++ )
  102.     {
  103.     j = bitrv( k, m );
  104.     if( j > k )
  105.         { /* executed approx. n/2 times */
  106.         p = &x[2*k];
  107.         tr = *p++;
  108.         ti = *p;
  109.         q = &x[2*j+1];
  110.         *p = *q;
  111.         *(--p) = *(--q);
  112.         *q++ = tr;
  113.         *q = ti;
  114.         }
  115.     }
  116.     
  117. /*                            fftr.c    */
  118. /*            Radix 2 Complex FFT            */
  119. n2 = n/2;
  120. nj = 1;
  121. pth = 1;
  122. dth = 0;
  123. th = 0;
  124.  
  125. for( l=0; l<m; l++ )
  126.     {    /* executed log2(n) times, total */
  127.     j = 0;
  128.     do
  129.         {    /* executed n-1 times, total */
  130.         r = th << 1;
  131.         si = sine[r];
  132.         co = sine[ n4 - r ];
  133.         if( j >= pth )
  134.             {
  135.             th -= dth;
  136.             co = -co;
  137.             }
  138.         else
  139.             th += dth;
  140.  
  141.         nd = j;
  142.  
  143.         do
  144.             { /* executed n/2 log2(n) times, total */
  145.             r = (nd << 1) + (nj << 1);
  146.             p = &x[ r ];
  147.             xr = *p++;
  148.             xi = *p;
  149.             tr = xr * co + xi * si;
  150.             ti = xi * co - xr * si;
  151.             r = nd << 1;
  152.             q = &x[ r ];
  153.             xr = *q++;
  154.             xi = *q;
  155.             *p = xi - ti;
  156.             *(--p) = xr - tr;
  157.             *q = xi + ti;
  158.             *(--q) = xr + tr;
  159.             nd += nj << 1;
  160.             }
  161.         while( nd < n );
  162.         }
  163.     while( ++j < nj );
  164.  
  165.     n2 >>= 1;
  166.     dth = n2;
  167.     pth = nj;
  168.     nj <<= 1;
  169.     }
  170.  
  171. /*                            fftr.c    */
  172.  
  173. /*    Special trick algorithm            */
  174. /*    converts to spectrum of real series    */
  175.  
  176. /* Highest frequency term; add space to input array if wanted */
  177. /*
  178. x[2*n] = x[0] - x[1];
  179. x[2*n+1] = 0.0;
  180. */
  181.  
  182. /* Zero frequency term */
  183. x[0] = x[0] + x[1];
  184. x[1] = 0.0;
  185. n2 = n/2;
  186.  
  187. for( j=1; j<=n2; j++ )
  188.     {    /* executed n/2 times */
  189.     si = sine[j];
  190.     co = sine[ n4 - j ];
  191.     p = &x[ 2*j ];
  192.     xr = *p++;
  193.     xi = *p;
  194.     q = &x[ 2*(n-j) ];
  195.     tr = *q++;
  196.     ti = *q;
  197.     a = xr + tr;
  198.     b = xi + ti;
  199.     c = xr - tr;
  200.     d = xi - ti;
  201.     bc = b * co;
  202.     cs = c * si;
  203.     bs = b * si;
  204.     cc = c * co;
  205.     *p = ( d - bs - cc )/2.0;
  206.     *(--p) = ( a + bc - cs )/2.0;
  207.     *q = -( d + bs + cc )/2.0;
  208.     *(--q) = ( a - bc + cs )/2.0;
  209.     }
  210.  
  211. return(0);
  212. }
  213.  
  214. /*                            fftr.c    */
  215.  
  216. /*    Bit reverser    */
  217.  
  218. int bitrv( j, m )
  219. int j, m;
  220. {
  221. register int j1, ans;
  222. short k;
  223.  
  224. ans = 0;
  225. j1 = j;
  226.  
  227. for( k=0; k<m; k++ )
  228.     {
  229.     ans = (ans << 1) + (j1 & 1);
  230.     j1 >>= 1;
  231.     }
  232.  
  233. return( ans );
  234. }
  235.